{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7cb97169",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyodeint\n",
    "import joblib\n",
    "import pandas as pd\n",
    "import random\n",
    "from chempy import Substance, Reaction, ReactionSystem\n",
    "from chempy.kinetics.rates import Arrhenius, MassAction\n",
    "from collections import defaultdict\n",
    "from chempy.units import SI_base_registry, default_units as u\n",
    "import numpy as np\n",
    "from chempy.kinetics.ode import get_odesys\n",
    "from itertools import product\n",
    "import csv\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "bf3e5ed6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Number of Reactions</th>\n",
       "      <th>Step</th>\n",
       "      <th>SM_left</th>\n",
       "      <th>C_left</th>\n",
       "      <th>B_left</th>\n",
       "      <th>P_right</th>\n",
       "      <th>IMP1_right</th>\n",
       "      <th>INT1_right</th>\n",
       "      <th>C_right</th>\n",
       "      <th>B_right</th>\n",
       "      <th>R_left</th>\n",
       "      <th>P_left</th>\n",
       "      <th>INT1_left</th>\n",
       "      <th>IMP2_right</th>\n",
       "      <th>Reaction Format</th>\n",
       "      <th>Initial Compounds</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>SM -&gt; INT1</td>\n",
       "      <td>SM</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>SM + R -&gt; INT1</td>\n",
       "      <td>SM,R</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>SM -&gt; IMP1</td>\n",
       "      <td>SM,IMP1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>SM + R -&gt; IMP1</td>\n",
       "      <td>SM,R,IMP1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>SM -&gt; IMP1 + INT1</td>\n",
       "      <td>SM,IMP1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>585</th>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>SM + P + C + B + R + INT1 -&gt; IMP2 + B</td>\n",
       "      <td>SM,B,C,R</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>586</th>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>SM + P + C + B + INT1 -&gt; IMP2 + C</td>\n",
       "      <td>SM,B,C</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>587</th>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>SM + P + C + B + R + INT1 -&gt; IMP2 + C</td>\n",
       "      <td>SM,B,C,R</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>588</th>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>SM + P + C + B + INT1 -&gt; IMP2 + C + B</td>\n",
       "      <td>SM,B,C</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>589</th>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>SM + P + C + B + R + INT1 -&gt; IMP2 + C + B</td>\n",
       "      <td>SM,B,C,R</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>590 rows × 16 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "     Number of Reactions  Step  SM_left  C_left  B_left  P_right  IMP1_right  \\\n",
       "0                      2     1        1       0       0        0           0   \n",
       "1                      2     1        1       0       0        0           0   \n",
       "2                      2     1        1       0       0        0           1   \n",
       "3                      2     1        1       0       0        0           1   \n",
       "4                      2     1        1       0       0        0           1   \n",
       "..                   ...   ...      ...     ...     ...      ...         ...   \n",
       "585                    2     2        1       1       1        0           0   \n",
       "586                    2     2        1       1       1        0           0   \n",
       "587                    2     2        1       1       1        0           0   \n",
       "588                    2     2        1       1       1        0           0   \n",
       "589                    2     2        1       1       1        0           0   \n",
       "\n",
       "     INT1_right  C_right  B_right  R_left  P_left  INT1_left  IMP2_right  \\\n",
       "0             1        0        0       0       0          0           0   \n",
       "1             1        0        0       1       0          0           0   \n",
       "2             0        0        0       0       0          0           0   \n",
       "3             0        0        0       1       0          0           0   \n",
       "4             1        0        0       0       0          0           0   \n",
       "..          ...      ...      ...     ...     ...        ...         ...   \n",
       "585           0        0        1       1       1          1           1   \n",
       "586           0        1        0       0       1          1           1   \n",
       "587           0        1        0       1       1          1           1   \n",
       "588           0        1        1       0       1          1           1   \n",
       "589           0        1        1       1       1          1           1   \n",
       "\n",
       "                               Reaction Format Initial Compounds  \n",
       "0                                   SM -> INT1                SM  \n",
       "1                               SM + R -> INT1              SM,R  \n",
       "2                                   SM -> IMP1           SM,IMP1  \n",
       "3                               SM + R -> IMP1         SM,R,IMP1  \n",
       "4                            SM -> IMP1 + INT1           SM,IMP1  \n",
       "..                                         ...               ...  \n",
       "585      SM + P + C + B + R + INT1 -> IMP2 + B          SM,B,C,R  \n",
       "586          SM + P + C + B + INT1 -> IMP2 + C            SM,B,C  \n",
       "587      SM + P + C + B + R + INT1 -> IMP2 + C          SM,B,C,R  \n",
       "588      SM + P + C + B + INT1 -> IMP2 + C + B            SM,B,C  \n",
       "589  SM + P + C + B + R + INT1 -> IMP2 + C + B          SM,B,C,R  \n",
       "\n",
       "[590 rows x 16 columns]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "database = joblib.load(\"database_two_reactions.joblib\")\n",
    "database"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "41252f1a",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(0, 146), (0, 147), (0, 154), (0, 155), (1, 146), (1, 147), (1, 154), (1, 155), (4, 146), (4, 147), (4, 154), (4, 155), (5, 146), (5, 147), (5, 154), (5, 155), (7, 126), (7, 162), (7, 163), (7, 334), (7, 382), (8, 138), (8, 170), (9, 138), (9, 139), (9, 170), (9, 171), (11, 126), (11, 162), (11, 163), (11, 334), (11, 382), (12, 138), (12, 170), (13, 138), (13, 139), (13, 170), (13, 171), (16, 146), (16, 147), (16, 154), (16, 155), (16, 200), (16, 201), (16, 208), (16, 209), (17, 146), (17, 147), (17, 154), (17, 155), (17, 200), (17, 201), (17, 208), (17, 209), (24, 146), (24, 147), (24, 154), (24, 155), (24, 200), (24, 201), (24, 208), (24, 209), (25, 146), (25, 147), (25, 154), (25, 155), (25, 200), (25, 201), (25, 208), (25, 209), (29, 126), (29, 162), (29, 163), (29, 179), (29, 216), (29, 217), (29, 334), (29, 382), (29, 400), (29, 448), (32, 138), (32, 170), (32, 192), (32, 224), (33, 138), (33, 139), (33, 170), (33, 171), (33, 192), (33, 193), (33, 224), (33, 225), (37, 126), (37, 162), (37, 163), (37, 179), (37, 216), (37, 217), (37, 334), (37, 382), (37, 400), (37, 448), (40, 138), (40, 170), (40, 192), (40, 224), (41, 138), (41, 139), (41, 170), (41, 171), (41, 192), (41, 193), (41, 224), (41, 225), (42, 150), (42, 151), (42, 158), (42, 159), (43, 150), (43, 151), (43, 158), (43, 159), (44, 146), (44, 147), (44, 154), (44, 155), (44, 254), (44, 255), (44, 262), (44, 263), (45, 146), (45, 147), (45, 154), (45, 155), (45, 254), (45, 255), (45, 262), (45, 263), (50, 150), (50, 151), (50, 158), (50, 159), (51, 150), (51, 151), (51, 158), (51, 159), (52, 146), (52, 147), (52, 154), (52, 155), (52, 254), (52, 255), (52, 262), (52, 263), (53, 146), (53, 147), (53, 154), (53, 155), (53, 254), (53, 255), (53, 262), (53, 263), (55, 128), (55, 166), (55, 167), (55, 338), (55, 386), (57, 126), (57, 162), (57, 163), (57, 232), (57, 270), (57, 271), (57, 334), (57, 382), (57, 466), (57, 514), (58, 142), (58, 174), (59, 142), (59, 143), (59, 174), (59, 175), (60, 138), (60, 170), (60, 246), (60, 278), (61, 138), (61, 139), (61, 170), (61, 171), (61, 246), (61, 247), (61, 278), (61, 279), (63, 128), (63, 166), (63, 167), (63, 338), (63, 386), (65, 126), (65, 162), (65, 163), (65, 232), (65, 270), (65, 271), (65, 334), (65, 382), (65, 466), (65, 514), (66, 142), (66, 174), (67, 142), (67, 143), (67, 174), (67, 175), (68, 138), (68, 170), (68, 246), (68, 278), (69, 138), (69, 139), (69, 170), (69, 171), (69, 246), (69, 247), (69, 278), (69, 279), (72, 150), (72, 151), (72, 158), (72, 159), (72, 204), (72, 205), (72, 212), (72, 213), (73, 150), (73, 151), (73, 158), (73, 159), (73, 204), (73, 205), (73, 212), (73, 213), (76, 146), (76, 147), (76, 154), (76, 155), (76, 200), (76, 201), (76, 208), (76, 209), (76, 254), (76, 255), (76, 262), (76, 263), (76, 308), (76, 309), (76, 316), (76, 317), (77, 146), (77, 147), (77, 154), (77, 155), (77, 200), (77, 201), (77, 208), (77, 209), (77, 254), (77, 255), (77, 262), (77, 263), (77, 308), (77, 309), (77, 316), (77, 317), (88, 150), (88, 151), (88, 158), (88, 159), (88, 204), (88, 205), (88, 212), (88, 213), (89, 150), (89, 151), (89, 158), (89, 159), (89, 204), (89, 205), (89, 212), (89, 213), (92, 146), (92, 147), (92, 154), (92, 155), (92, 200), (92, 201), (92, 208), (92, 209), (92, 254), (92, 255), (92, 262), (92, 263), (92, 308), (92, 309), (92, 316), (92, 317), (93, 146), (93, 147), (93, 154), (93, 155), (93, 200), (93, 201), (93, 208), (93, 209), (93, 254), (93, 255), (93, 262), (93, 263), (93, 308), (93, 309), (93, 316), (93, 317), (97, 128), (97, 166), (97, 167), (97, 181), (97, 220), (97, 221), (97, 338), (97, 386), (97, 404), (97, 452), (101, 126), (101, 162), (101, 163), (101, 179), (101, 216), (101, 217), (101, 232), (101, 270), (101, 271), (101, 285), (101, 324), (101, 325), (101, 334), (101, 382), (101, 400), (101, 448), (101, 466), (101, 514), (101, 532), (101, 580), (104, 142), (104, 174), (104, 196), (104, 228), (105, 142), (105, 143), (105, 174), (105, 175), (105, 196), (105, 197), (105, 228), (105, 229), (108, 138), (108, 170), (108, 192), (108, 224), (108, 246), (108, 278), (108, 300), (108, 332), (109, 138), (109, 139), (109, 170), (109, 171), (109, 192), (109, 193), (109, 224), (109, 225), (109, 246), (109, 247), (109, 278), (109, 279), (109, 300), (109, 301), (109, 332), (109, 333), (113, 128), (113, 166), (113, 167), (113, 181), (113, 220), (113, 221), (113, 338), (113, 386), (113, 404), (113, 452), (117, 126), (117, 162), (117, 163), (117, 179), (117, 216), (117, 217), (117, 232), (117, 270), (117, 271), (117, 285), (117, 324), (117, 325), (117, 334), (117, 382), (117, 400), (117, 448), (117, 466), (117, 514), (117, 532), (117, 580), (120, 142), (120, 174), (120, 196), (120, 228), (121, 142), (121, 143), (121, 174), (121, 175), (121, 196), (121, 197), (121, 228), (121, 229), (124, 138), (124, 170), (124, 192), (124, 224), (124, 246), (124, 278), (124, 300), (124, 332), (125, 138), (125, 139), (125, 170), (125, 171), (125, 192), (125, 193), (125, 224), (125, 225), (125, 246), (125, 247), (125, 278), (125, 279), (125, 300), (125, 301), (125, 332), (125, 333)]\n",
      "456\n"
     ]
    }
   ],
   "source": [
    "valid_combinations_two_reactions = joblib.load(\"valid_combinations_two_reactions.joblib\")\n",
    "print(valid_combinations_two_reactions)\n",
    "print(len(valid_combinations_two_reactions))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "f6fb6553",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1, 2, 3]\n"
     ]
    }
   ],
   "source": [
    "equiv_B = [i + 1 for i in range(3)]\n",
    "print(equiv_B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "ca9e7fe6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "8\n"
     ]
    }
   ],
   "source": [
    "# Generate initial conditions\n",
    "seed_value = 37\n",
    "random.seed(seed_value)\n",
    "\n",
    "T_list = [273, 313, 373]\n",
    "\n",
    "conc_SM = [0.3] \n",
    "conc_IMP2 = 0\n",
    "conc_IMP3 = 0\n",
    "conc_IMP4 = 0\n",
    "\n",
    "equiv_R = [1.2, 1.5]\n",
    "equiv_B = [3] \n",
    "equiv_C = [0.3] \n",
    "equiv_P = [0, 0.5]\n",
    "equiv_IMP1 = [0, 0.002]\n",
    "\n",
    "imp2_array = np.array(conc_IMP2)\n",
    "imp3_array = np.array(conc_IMP3)\n",
    "imp4_array = np.array(conc_IMP4)\n",
    "\n",
    "c0_list = []\n",
    "for conc_sm in conc_SM:\n",
    "    for equiv_r in equiv_R:\n",
    "        for equiv_b in equiv_B:\n",
    "            for equiv_c in equiv_C:\n",
    "                for equiv_p in equiv_P:\n",
    "                    for equiv_imp1 in equiv_IMP1:\n",
    "                        sm_array = np.array(conc_sm)\n",
    "                        r_array = np.array(equiv_r*conc_sm)\n",
    "                        b_array = np.array(equiv_b*conc_sm)\n",
    "                        c_array = np.array(equiv_c*conc_sm)\n",
    "                        p_array = np.array(equiv_p*conc_sm)\n",
    "                        imp1_array = np.array(equiv_imp1*conc_sm)\n",
    "                        c0_list.append(\"c0 = defaultdict(lambda: 0*u.M, {{'SM': {}*u.M, 'R': {}*u.M, 'B': {}*u.M, 'C': {}*u.M, 'P': {}*u.M, 'IMP1': {}*u.M, 'IMP2': {}*u.M}})\"\\\n",
    "                                       .format(sm_array, r_array, b_array, c_array, p_array, imp1_array, imp2_array))\n",
    "print(len(c0_list))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "f997d27c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[\"c0 = defaultdict(lambda: 0*u.M, {'SM': 0.3*u.M, 'R': 0.36*u.M, 'B': 0.8999999999999999*u.M, 'C': 0.09*u.M, 'P': 0.0*u.M, 'IMP1': 0.0*u.M, 'IMP2': 0*u.M})\",\n",
       " \"c0 = defaultdict(lambda: 0*u.M, {'SM': 0.3*u.M, 'R': 0.36*u.M, 'B': 0.8999999999999999*u.M, 'C': 0.09*u.M, 'P': 0.0*u.M, 'IMP1': 0.0006*u.M, 'IMP2': 0*u.M})\",\n",
       " \"c0 = defaultdict(lambda: 0*u.M, {'SM': 0.3*u.M, 'R': 0.36*u.M, 'B': 0.8999999999999999*u.M, 'C': 0.09*u.M, 'P': 0.15*u.M, 'IMP1': 0.0*u.M, 'IMP2': 0*u.M})\",\n",
       " \"c0 = defaultdict(lambda: 0*u.M, {'SM': 0.3*u.M, 'R': 0.36*u.M, 'B': 0.8999999999999999*u.M, 'C': 0.09*u.M, 'P': 0.15*u.M, 'IMP1': 0.0006*u.M, 'IMP2': 0*u.M})\",\n",
       " \"c0 = defaultdict(lambda: 0*u.M, {'SM': 0.3*u.M, 'R': 0.44999999999999996*u.M, 'B': 0.8999999999999999*u.M, 'C': 0.09*u.M, 'P': 0.0*u.M, 'IMP1': 0.0*u.M, 'IMP2': 0*u.M})\",\n",
       " \"c0 = defaultdict(lambda: 0*u.M, {'SM': 0.3*u.M, 'R': 0.44999999999999996*u.M, 'B': 0.8999999999999999*u.M, 'C': 0.09*u.M, 'P': 0.0*u.M, 'IMP1': 0.0006*u.M, 'IMP2': 0*u.M})\",\n",
       " \"c0 = defaultdict(lambda: 0*u.M, {'SM': 0.3*u.M, 'R': 0.44999999999999996*u.M, 'B': 0.8999999999999999*u.M, 'C': 0.09*u.M, 'P': 0.15*u.M, 'IMP1': 0.0*u.M, 'IMP2': 0*u.M})\",\n",
       " \"c0 = defaultdict(lambda: 0*u.M, {'SM': 0.3*u.M, 'R': 0.44999999999999996*u.M, 'B': 0.8999999999999999*u.M, 'C': 0.09*u.M, 'P': 0.15*u.M, 'IMP1': 0.0006*u.M, 'IMP2': 0*u.M})\"]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c0_list"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "851aa765",
   "metadata": {},
   "source": [
    "### Functions Used for Generating Time-Course Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "1a6f36be",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define a function to convert a flat combination into dictionaries\n",
    "def flat_to_dicts(flat_combination, reagents, products):\n",
    "    reagent_values = flat_combination[:len(reagents)]\n",
    "    product_values = flat_combination[len(reagents):]\n",
    "    \n",
    "    reagents_dict = {reagent: value for reagent, value in zip(reagents, reagent_values)}\n",
    "    products_dict = {product: value for product, value in zip(products, product_values)}\n",
    "    \n",
    "    return reagents_dict, products_dict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "bdb806f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert list to daframe: Initialize dataframe and add data to it\n",
    "def generate_df_time_course_data(data):\n",
    "   con_SM, con_R, con_B, con_C, con_P, con_IMP1, con_IMP2, con_INT1 = \\\n",
    "      [], [], [], [], [], [], [], []\n",
    "   timestamp = [0,0.1,1,20,40,60,120,180,240,300,360,480,600,720,840,960,1080,1200,\\\n",
    "                1500,1800,2400,3000,3600,4500,5400,6300,7200,9000,10800,14400]\n",
    "   for i in range(0, len(timestamp)):\n",
    "      con_SM.append(\"cSM_\" + str(timestamp[i]) + \"s\") \n",
    "      con_R.append(\"cR_\" + str(timestamp[i]) + \"s\") \n",
    "      con_B.append(\"cB_\" + str(timestamp[i]) + \"s\")\n",
    "      con_C.append(\"cC_\" + str(timestamp[i]) + \"s\")\n",
    "      con_P.append(\"cP_\" + str(timestamp[i]) + \"s\") \n",
    "      con_IMP1.append(\"cIMP1_\" + str(timestamp[i]) + \"s\") \n",
    "      con_IMP2.append(\"cIMP2_\" + str(timestamp[i]) + \"s\")\n",
    "      con_INT1.append(\"cINT1_\" + str(timestamp[i]) + \"s\") \n",
    "      \n",
    "   columns = ['Temperature', 'A1', 'Ea1', 'A2', 'Ea2', 'A3', 'Ea3', 'A4', 'Ea4'] + \\\n",
    "      con_SM + con_R + con_B + con_C + con_P + con_IMP1 + con_IMP2 +con_INT1 +  \\\n",
    "         ['Fast_rxn1', 'Medium_rxn1', 'Slow_rxn1', 'Fast_rxn2', 'Medium_rxn2', 'Slow_rxn2',  'Reaction_order', 'Mechanism']\n",
    "   df_time_course_data = pd.DataFrame(data, columns = columns)\n",
    "   return df_time_course_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "982541bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_order_combinations(reagents, products):\n",
    "    values_reagents = [1, 2]\n",
    "\n",
    "    # Generate all possible combinations of values for reagents. The reaction order of product is always 1.\n",
    "    combinations_reagents = product(values_reagents, repeat=len(reagents))\n",
    "    product_combination = tuple([1] * len(products))\n",
    "    \n",
    "    # Initialize an empty list to store the combinations\n",
    "    combinations = []\n",
    "    # Generate combinations for reagents\n",
    "    for reagent_combination in combinations_reagents:\n",
    "        # Since products should always be 1, create a tuple of 1s for products\n",
    "        product_combination = tuple([1] * len(products))\n",
    "        \n",
    "        # Combine reagent and product combinations\n",
    "        combination = list(reagent_combination) + list(product_combination)\n",
    "        combinations.append(combination)\n",
    "    return(combinations)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "7c86b38a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def initialize_csv(filename):\n",
    "    with open('csv_header.csv', 'r') as header_file:\n",
    "        csvreader = csv.reader(header_file)\n",
    "        header = next(csvreader) \n",
    "    with open(filename, 'w', newline='') as csvfile:\n",
    "        csvwriter = csv.writer(csvfile)\n",
    "        csvwriter.writerow(header)  \n",
    "\n",
    "def write_row_to_csv(filename, row):\n",
    "    with open(filename, 'a', newline='') as csvfile:\n",
    "        csvwriter = csv.writer(csvfile)\n",
    "        csvwriter.writerow(row)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f382cc2",
   "metadata": {},
   "source": [
    "### Generate Time Course Data for Two-reaction Mechnism"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "43781ae3",
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_time_course_data_two_rxn_csv_out_sim(two_reaction, c0, T, csv_filename, max_file_size=50*1024*1024):\n",
    "\n",
    "    if not os.path.isfile(csv_filename):\n",
    "        initialize_csv(csv_filename)\n",
    "    \n",
    "    local_vars = {}\n",
    "    exec(c0, globals(), local_vars)\n",
    "    c0 = local_vars['c0']\n",
    "\n",
    "    counter = 0\n",
    "\n",
    "    for mechanism in two_reaction:\n",
    "        original_list = database.iloc[mechanism[0],-2].split(' -> ')\n",
    "        reagents_rxn1 = original_list[0].split(' + ')\n",
    "        products_rxn1 = original_list[1].split(' + ')\n",
    "        original_list = database.iloc[mechanism[1],-2].split(' -> ')\n",
    "        reagents_rxn2 = original_list[0].split(' + ')\n",
    "        products_rxn2 = original_list[1].split(' + ')\n",
    "        # print(reagents_rxn1)\n",
    "        # print(products_rxn1)\n",
    "        # print(reagents_rxn2)\n",
    "        # print(products_rxn2)\n",
    "        \n",
    "        # Generate all possible combinations for the reaction order\n",
    "        # Create a list of all possible values for reaction order (1 and 2)\n",
    "        values = [1, 2]\n",
    "\n",
    "        # Generate all possible combinations of values for reagents. The reaction order of product is always 1.\n",
    "        combinations_rxn1 = generate_order_combinations(reagents_rxn1, products_rxn1)\n",
    "        combinations_rxn2 = generate_order_combinations(reagents_rxn2, products_rxn2)\n",
    "\n",
    "    \n",
    "        # Find the order of species in time course data\n",
    "        all_species = reagents_rxn1 + products_rxn1 + reagents_rxn2 + products_rxn2\n",
    "        # Merge all lists and remove duplicates\n",
    "        merged_list = list(set(all_species))\n",
    "        # Sort the merged list alphabetically\n",
    "        merged_list.sort()\n",
    "        # Create a dictionary with items as keys and their corresponding indexes as values\n",
    "        species_order = {item: index for index, item in enumerate(merged_list)}\n",
    "\n",
    "        A1, A2, A3, A4 = [0,0,0,0]\n",
    "        Ea1, Ea2, Ea3, Ea4 = [0,0,0,0]\n",
    "        timestamp = [0,0.1,1,20,40,60,120,180,240,300,360,480,600,720,840,960,1080,1200,1500,1800,2400,3000,3600,4500,5400,6300,7200,9000,10800,14400]\n",
    "        num_rxn = 3\n",
    "        \n",
    "        for i, combination1 in enumerate(combinations_rxn1):\n",
    "            for j, combination2 in enumerate(combinations_rxn2):\n",
    "                # Generate randome Ea and A1\n",
    "                # 10 time course data for each category (low, medium, fast)\n",
    "                for a in range(0, num_rxn):\n",
    "                    for b in range(0, num_rxn):\n",
    "                        Ea1 = random.randint(0,200)\n",
    "                        Ea2 = random.randint(0,200)\n",
    "                        rate = [0] * 6\n",
    "                        # For slow reactions\n",
    "                        if a < num_rxn/3:\n",
    "                            A1 = random.uniform(10e-5, 0.1)\n",
    "                            rate[2] = 1\n",
    "                        if a >= num_rxn/3 and a < 2*num_rxn/3:\n",
    "                            A1 = random.uniform(0.1, 10)\n",
    "                            rate[1] = 1\n",
    "                        if a >= 2*num_rxn/3:\n",
    "                            A1 = random.uniform(10, 1000)\n",
    "                            rate[0] = 1\n",
    "                        if b < num_rxn/3:\n",
    "                            A2 = random.uniform(10e-5, 0.1)\n",
    "                            rate[5] = 1\n",
    "                        if b >= num_rxn/3 and b < 2*num_rxn/3:\n",
    "                            A2 = random.uniform(0.1, 10)\n",
    "                            rate[4] = 1\n",
    "                        if b >= 2*num_rxn/3:\n",
    "                            A2 = random.uniform(10, 1000)\n",
    "                            rate[3] = 1\n",
    "\n",
    "                        reagents_dict1, products_dict1 = flat_to_dicts(combination1, reagents_rxn1, products_rxn1)\n",
    "                        reaction1 = f\"r{b}_1 = Reaction({reagents_dict1}, {products_dict1}, MassAction(Arrhenius(unique_keys=('A1', 'Ea_R_1'))))\"\n",
    "                        # print(reaction1)\n",
    "                        exec(reaction1)\n",
    "                        \n",
    "                        reagents_dict2, products_dict2 = flat_to_dicts(combination2, reagents_rxn2, products_rxn2)\n",
    "                        reaction2 = f\"r{b}_2 = Reaction({reagents_dict2}, {products_dict2}, MassAction(Arrhenius(unique_keys=('A2', 'Ea_R_2'))))\"\n",
    "                        # print(reaction2)\n",
    "                        exec(reaction2)\n",
    "\n",
    "                        mech = [eval(\"r%s_1\" % b), eval(\"r%s_2\" % b)]\n",
    "                        rsys = ReactionSystem(mech)\n",
    "\n",
    "                        # Generate time course data\n",
    "                        variables = c0.copy()\n",
    "                        order_rxn1 = eval(\"r%s_1.order()\" % b)\n",
    "                        order_rxn2 = eval(\"r%s_2.order()\" % b)\n",
    "                        unit_orderofmag_rxn1 = order_rxn1 - 1\n",
    "                        unit_orderofmag_rxn2 = order_rxn2 - 1\n",
    "                        params = {'A1': A1/u.second*u.M**-unit_orderofmag_rxn1, 'Ea_R_1': Ea1/8.3145*u.K**-1, \\\n",
    "                                'A2': A2/u.second*u.M**-unit_orderofmag_rxn2, 'Ea_R_2': Ea2/8.3145*u.K**-1, 'temperature': T*u.K}\n",
    "\n",
    "                        variables.update(params)\n",
    "                        rsys.rates(variables)\n",
    "                        odesys, extra = get_odesys(rsys, include_params=False, lower_bounds=0)\n",
    "                        # results = odesys.integrate(4*3600*u.s, c0, params, integrator='odeint')\n",
    "                        try:\n",
    "                            results = odesys.integrate(4 * 3600 * u.s, c0, params, integrator='odeint', atol=1e-6, rtol=1e-6, nsteps = 5000)\n",
    "                            conc = results.at(timestamp)\n",
    "                            res = [T, A1, Ea1, A2, Ea2, A3, Ea3, A4, Ea4]\n",
    "                            \n",
    "                            # Append concnetrations\n",
    "                            all_species = ['SM', 'R', 'B', 'C', 'P', 'IMP1', 'IMP2', 'INT1']\n",
    "                            for species in all_species:\n",
    "                                if species in species_order:\n",
    "                                    for m in range(0, len(timestamp)):\n",
    "                                        # Remove negative values\n",
    "                                        if conc[m][0][species_order[species]] > 0:\n",
    "                                            res.append(conc[m][0][species_order[species]])\n",
    "                                        else:\n",
    "                                            res.append(0)\n",
    "                                else:\n",
    "                                    res += ([0]* len(timestamp))\n",
    "\n",
    "                            rxn_order = {'rxn1': reagents_dict1, 'rxn_2': reagents_dict2}\n",
    "                            res += rate + [rxn_order] + ([mechanism])\n",
    "\n",
    "                            # write the csv\n",
    "                            counter += 1\n",
    "                            write_row_to_csv(csv_filename, res)\n",
    "                            file_size = os.path.getsize(csv_filename)\n",
    "                            # print(file_size)\n",
    "                            # print(max_file_size)\n",
    "                            if file_size >= max_file_size:    \n",
    "                                print('yes')\n",
    "                                import datetime\n",
    "                                current_datetime = datetime.datetime.now()\n",
    "                                formatted_datetime = current_datetime.strftime(\"%Y-%m-%d_%H%M%S\")\n",
    "                                new_filename = f\"two_reactions_{formatted_datetime}.csv\"\n",
    "                                # counter +=1\n",
    "                                initialize_csv(new_filename)\n",
    "                                csv_filename = new_filename\n",
    "                        except RuntimeError:\n",
    "                            print('pass')\n",
    "                            pass\n",
    "                        \n",
    "    print('Number of time course data generated: ', counter)\n",
    "    # return(data_two_rxn)\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "8ec0be56",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "Number of time course data generated:  240408\n"
     ]
    }
   ],
   "source": [
    "generate_time_course_data_two_rxn_csv_out_sim(valid_combinations_two_reactions, c0_list[0], T_list[0], 'two_reactions_273K_1.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "a17d38a3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "Number of time course data generated:  240408\n"
     ]
    }
   ],
   "source": [
    "generate_time_course_data_two_rxn_csv_out_sim(valid_combinations_two_reactions, c0_list[1], T_list[0], 'two_reactions_273K_2.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "2f4ee9e4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "Number of time course data generated:  240408\n"
     ]
    }
   ],
   "source": [
    "generate_time_course_data_two_rxn_csv_out_sim(valid_combinations_two_reactions, c0_list[2], T_list[0], 'two_reactions_273K_3.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "506ad66b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "yes\n",
      "yes\n",
      "yes\n",
      "pass\n",
      "pass\n",
      "yes\n",
      "pass\n",
      "pass\n",
      "pass\n",
      "pass\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "pass\n",
      "pass\n",
      "pass\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "pass\n",
      "pass\n",
      "pass\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "pass\n",
      "pass\n",
      "pass\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "Number of time course data generated:  240393\n"
     ]
    }
   ],
   "source": [
    "generate_time_course_data_two_rxn_csv_out_sim(valid_combinations_two_reactions, c0_list[3], T_list[0], 'two_reactions_273K_4.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "ec660ff0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "yes\n",
      "yes\n",
      "yes\n",
      "pass\n",
      "pass\n",
      "yes\n",
      "pass\n",
      "pass\n",
      "pass\n",
      "pass\n",
      "pass\n",
      "yes\n",
      "yes\n",
      "pass\n",
      "pass\n",
      "pass\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "pass\n",
      "pass\n",
      "pass\n",
      "pass\n",
      "pass\n",
      "yes\n",
      "pass\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "pass\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "Number of time course data generated:  240391\n"
     ]
    }
   ],
   "source": [
    "generate_time_course_data_two_rxn_csv_out_sim(valid_combinations_two_reactions, c0_list[4], T_list[0], 'two_reactions_273K_5.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "5b36eb76",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "yes\n",
      "yes\n",
      "pass\n",
      "pass\n",
      "pass\n",
      "yes\n",
      "pass\n",
      "pass\n",
      "pass\n",
      "yes\n",
      "pass\n",
      "pass\n",
      "pass\n",
      "yes\n",
      "yes\n",
      "pass\n",
      "pass\n",
      "pass\n",
      "pass\n",
      "pass\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "pass\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "pass\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "Number of time course data generated:  240392\n"
     ]
    }
   ],
   "source": [
    "generate_time_course_data_two_rxn_csv_out_sim(valid_combinations_two_reactions, c0_list[5], T_list[0], 'two_reactions_273K_6.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "345b44b7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pass\n",
      "pass\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "Number of time course data generated:  240406\n"
     ]
    }
   ],
   "source": [
    "generate_time_course_data_two_rxn_csv_out_sim(valid_combinations_two_reactions, c0_list[6], T_list[0], 'two_reactions_273K_7.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "e8f06136",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "yes\n",
      "pass\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "yes\n",
      "Number of time course data generated:  240407\n"
     ]
    }
   ],
   "source": [
    "generate_time_course_data_two_rxn_csv_out_sim(valid_combinations_two_reactions, c0_list[7], T_list[0], 'two_reactions_273K_8.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "507044ee",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
